# setwd("~/Dropbox/T&G project survey/Ukraine/replication_JOP")
# libraries
library(tidyverse)
library(haven)
library(cobalt)
library(patchwork)
library(paletteer)
library(MatchIt)
library(estimatr)
library(lmtest)
library(stargazer)
library(xtable)
source("func.R")

# ggplot defaults: geoms theme
update_geom_defaults("text", list(family = "Archivo Narrow"))
update_geom_defaults("label", list(family = "Archivo Narrow"))
theme_set(theme_nice())

## Load data and transform
d = read_dta("data/bbdd_voto_fechas.dta", encoding = "latin1") %>%
  mutate(
    pr_voto = ifelse(p1 == 99, NA, p1),
    pr_voto_psoe = ifelse(p4_1 == 99, NA, p4_1),
    pr_voto_pp = ifelse(p4_2 == 99, NA, p4_2),
    pr_voto_up = ifelse(p4_3 == 99, NA, p4_3),
    pr_voto_vox = ifelse(p4_4 == 99, NA, p4_4),
    ideology = ifelse(p6 == 99, NA, p6)) %>%
  mutate(post = ifelse(fecha >= 24, 1, 0)) %>%
  mutate(date = as.Date(paste0("2022-02-", fecha))) %>%
  mutate(
    voto_psoe = ifelse(p3 == 1, 1, 0),
    voto_pp = ifelse(p3 == 2, 1, 0),
    voto_up = ifelse(p3 == 4, 1, 0),
    voto_vox = ifelse(p3 == 5, 1, 0)) %>%
  filter(!is.na(ideology))

# Dictionary
dict = tribble(~dirty, ~clean,
  "sexo_2", "Sex (female)",
  "edad", "Age",
  "clase_social_r", "Social class\n(scale)",
  "hab_r", "Location\n(below 10,000 residents)",
  "educacion_r", "Education\n(scale)",
  "Q12_2", "Voted in last elections",
  "Q14", "Income\n(scale)",
  "ideology", "Ideology\n(scale)")

# Dates
dates = ggplot(d, aes(x = date)) +
  geom_histogram(binwidth = 1, position = "dodge2", color = "white") +
  theme(legend.position = "top") +
  scale_x_date(date_breaks = "1 day", date_labels = "%b %d") +
  scale_y_continuous(limits = c(0, 390)) +
  geom_segment(
    aes(x = as.Date("2022-02-24")-0.25, y = 0, xend = as.Date("2022-02-24")-0.25, yend = 350),
    linetype = "dashed", color = "red") +
  annotate(geom = "text", x = as.Date("2022-02-24")-0.15, color = "red", hjust = 0,
    y = 350, label = "February 24 (~5h CEST):\nRussian invasion starts") +
  labs(x = "", y = "Interviews per day")
ggsave("figures/hist_dates_barometer.pdf", height = 3, width = 5, device = cairo_pdf)

## Check balance

pdat_diff = love.plot(post ~ sexo + edad + clase_social_r + ideology +
    hab_r + educacion_r, data = d,
    thresholds = c(m = .1, v = 2))$data %>%
  tibble() %>%
  left_join(dict, by = c("var" = "dirty")) %>%
  mutate(clean = factor(clean)) %>%
  mutate(balanced = ifelse(abs(stat) < .1, "Balanced", "Unbalanced")) %>%
  mutate(statistic = "Standardized mean differences")

balance = pdat_diff %>%
  ggplot(aes(y = clean, x = stat, color = balanced)) +
  geom_point(size = 2.2) +
  geom_vline(xintercept = 0) +
  geom_vline(xintercept = -.1, lty = 2) +
  geom_vline(xintercept = .1, lty = 2) +
  coord_cartesian(xlim = c(-1, 1)) +
  labs(x = "Standardized mean differences", y = NULL,
       color = NULL) +
  theme(legend.position = "top") +
  scale_color_paletteer_d(`"wesanderson::Darjeeling1"`, direction = -1)

ggsave(balance, filename = "figures/barometer_pre_balance.pdf",
  width = 5, height = 5, device = cairo_pdf)


## Matching

# Two samples separately (with/without Feb 24)
f = as.formula("post ~ sexo + edad + clase_social_r + ideology + hab_r + educacion_r")
match_full = matchit(f, data = d, method = "full")
match_no24 = matchit(f, data = subset(d, fecha != 24), method = "full")
# Data
matched_data = match.data(match_full)
matched_data_no24 = match.data(match_no24)

# Balance table
bal_tab = rbind(
  matrix(rep(NA, 3), ncol=3),
  summary(match_full)$sum.all[, 1:3],
  matrix(rep(NA, 3), ncol=3),
  summary(match_full)$sum.matched[, 1:3]) %>%
  as.data.frame %>%
  mutate(var = gsub("\\.1", "", rownames(.))) %>%
  mutate(var = recode(var,
  "distance" = "Distance (propensity score)",
  "ideology" = "Ideology",
  "sexo" = "Sex",
  "edad" = "Age",
  "hab_r" = "Location\n(below 10,000 residents)",
  "clase_social_r" = "Social class",
  "educacion_r" = "Education")) %>%
  select(var, 1:3)

c = "Balance table for complementary survey"
l = "tab:balance_barometer"
bal_tab_tex = str_split(print(xtable(bal_tab, caption = c, label = l),
  include.rownames = FALSE, caption.placement = "top"), "\n")[[1]]
bal_tab_tex = gsub("X &  &  &  ", "\\\\textbf{Original dataset}&&&", bal_tab_tex)
bal_tab_tex[grepl("Original", bal_tab_tex)][2] = gsub("Original", "Matched",
  bal_tab_tex[grepl("Original", bal_tab_tex)][2])
bal_tab_tex = gsub("var ", "", bal_tab_tex)

filecon = file("tables/tab_balance_barometer.tex")
writeLines(bal_tab_tex, filecon)
close(filecon)

# Balance tables, excluding February 24th
bal_tab_no24 = rbind(
  matrix(rep(NA, 3), ncol=3),
  summary(match_no24)$sum.all[, 1:3],
  matrix(rep(NA, 3), ncol=3),
  summary(match_no24)$sum.matched[, 1:3]) %>%
  as.data.frame %>%
  mutate(var = gsub("\\.1", "", rownames(.))) %>%
  mutate(var = recode(var,
  "distance" = "Distance (propensity score)",
  "ideology" = "Ideology",
  "sexo" = "Sex",
  "edad" = "Age",
  "clase_social_r" = "Social class",
  "educacion_r" = "Education")) %>%
  select(var, 1:3)

c = "Balance table for complementary survey (excluding February 24)"
l = "tab:balance_barometer_no24"
bal_tab_tex_no24 = str_split(print(xtable(bal_tab_no24, caption = c, label = l),
  include.rownames = FALSE, caption.placement = "top"), "\n")[[1]]
bal_tab_tex_no24 = gsub("X &  &  &  ", "\\\\textbf{Original dataset}&&&", bal_tab_tex_no24)
bal_tab_tex_no24[grepl("Original", bal_tab_tex_no24)][2] = gsub("Original", "Matched",
  bal_tab_tex_no24[grepl("Original", bal_tab_tex_no24)][2])
bal_tab_tex_no24 = gsub("var ", "", bal_tab_tex_no24)

filecon = file("tables/tab_balance_barometer_no24.tex")
writeLines(bal_tab_tex_no24, filecon)
close(filecon)

# Run models and produce plots
m1 = lm_robust(pr_voto ~ post,
  data = matched_data, weights = weights, clusters = subclass, se_type = "stata")
m2 = lm_robust(pr_voto_psoe ~ post,
  data = matched_data, weights = weights, clusters = subclass, se_type = "stata")
m3 = lm_robust(pr_voto_pp ~ post,
  data = matched_data, weights = weights, clusters = subclass, se_type = "stata")
m4 = lm_robust(pr_voto_up ~ post,
  data = matched_data, weights = weights, clusters = subclass, se_type = "stata")
m5 = lm_robust(pr_voto_vox ~ post,
  data = matched_data, weights = weights, clusters = subclass, se_type = "stata")
m1b = lm_robust(pr_voto ~ post,
  data = matched_data_no24, weights = weights, clusters = subclass, se_type = "stata")
m2b = lm_robust(pr_voto_psoe ~ post,
  data = matched_data_no24, weights = weights, clusters = subclass, se_type = "stata")
m3b = lm_robust(pr_voto_pp ~ post,
  data = matched_data_no24, weights = weights, clusters = subclass, se_type = "stata")
m4b = lm_robust(pr_voto_up ~ post,
  data = matched_data_no24, weights = weights, clusters = subclass, se_type = "stata")
m5b = lm_robust(pr_voto_vox ~ post,
  data = matched_data_no24, weights = weights, clusters = subclass, se_type = "stata")

m_coefs = rbind(
    tidy(m1, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m1b, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2b, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3b, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4b, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5b, conf.int = TRUE) %>% filter(term == "post")) %>%
  mutate(outcome2 = recode(outcome,
    "pr_voto" = "Would vote",
    "pr_voto_psoe" = "Vote for PSOE",
    "pr_voto_pp" = "Vote for PP",
    "pr_voto_up" = "Vote for UP",
    "pr_voto_vox" = "Vote for VOX")) %>%
  mutate(model = rep(c("Keeping February 24th in sample",
    "Excluding February 24th from sample"), each = 5)) %>%
  mutate(model2 = recode(model, "Keeping February 24th in sample" = 1,
    "Excluding February 24th from sample" = 2)) %>%
  mutate(
    conf.high90 = estimate + std.error * qnorm(0.95),
    conf.low90 = estimate - std.error * qnorm(0.95))

main = ggplot(m_coefs, aes(x = estimate, y = outcome2)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
    fatten = 1, size = 2, shape = 21, fill = "white") +
  geom_pointrange(aes(xmin = conf.low90, xmax = conf.high90),
    fatten = 1, size = 3, shape = 21, fill = "white") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  labs(x = "\nCoefficient estimate of invasion", y = "",
    subtitle = "Bars indicate 90% and 95% CIs.") +
  theme(axis.text.y = element_text(size = 12)) +
  facet_wrap(~reorder(model, model2))
ggsave("figures/barometer_results.pdf", width = 10, height = 3, device = cairo_pdf)

## Regression tables, with/without controls

m1_t = lm(pr_voto ~ post, data = matched_data, weights = weights)
m2_t = lm(pr_voto_psoe ~ post, data = matched_data, weights = weights)
m3_t = lm(pr_voto_pp ~ post, data = matched_data, weights = weights)
m4_t = lm(pr_voto_up ~ post, data = matched_data, weights = weights)
m5_t = lm(pr_voto_vox ~ post, data = matched_data, weights = weights)
m1b_t = lm(pr_voto ~ post, data = matched_data_no24, weights = weights)
m2b_t = lm(pr_voto_psoe ~ post, data = matched_data_no24, weights = weights)
m3b_t = lm(pr_voto_pp ~ post, data = matched_data_no24, weights = weights)
m4b_t = lm(pr_voto_up ~ post, data = matched_data_no24, weights = weights)
m5b_t = lm(pr_voto_vox ~ post, data = matched_data_no24, weights = weights)

fright = "~ post + sexo + edad + clase_social_r + ideology + hab_r + educacion_r + factor(ccaa)"
m1_t_covs = lm(as.formula(paste0("pr_voto", fright)), data = d)
m2_t_covs = lm(as.formula(paste0("pr_voto_psoe", fright)), data = d)
m3_t_covs = lm(as.formula(paste0("pr_voto_pp", fright)), data = d)
m4_t_covs = lm(as.formula(paste0("pr_voto_up", fright)), data = d)
m5_t_covs = lm(as.formula(paste0("pr_voto_vox", fright)), data = d)
m1b_t_covs = lm(as.formula(paste0("pr_voto", fright)), data = subset(d, fecha != 24))
m2b_t_covs = lm(as.formula(paste0("pr_voto_psoe", fright)), data = subset(d, fecha != 24))
m3b_t_covs = lm(as.formula(paste0("pr_voto_pp", fright)), data = subset(d, fecha != 24))
m4b_t_covs = lm(as.formula(paste0("pr_voto_up", fright)), data = subset(d, fecha != 24))
m5b_t_covs = lm(as.formula(paste0("pr_voto_vox", fright)), data = subset(d, fecha != 24))

# Model lists
mlist = list(m1_t, m2_t, m3_t, m4_t, m5_t)
mlist_covs = list(m1_t_covs, m2_t_covs, m3_t_covs, m4_t_covs, m5_t_covs)
mlist_no24 = list(m1b_t, m2b_t, m3b_t, m4b_t, m5b_t)
mlist_no24_covs = list(m1b_t_covs, m2b_t_covs, m3b_t_covs, m4b_t_covs, m5b_t_covs)

# Dep var names
dvs = c("Would vote", "Vote for PSOE", "Vote for PP", "Vote for UP", "Vote for VOX")

# Table (full)
filecon = file("tables/tab_barometer_full.tex")
writeLines(
  stargazer(mlist,
    se = starprep(list(m1, m2, m3, m4, m5),
      clusters = matched_data$subclass, se_type = "stata"),
    keep = c("Constant", "post"),
    label = "tab:lm_baro",
    title = "Effect of invasion on electoral behavior",
    order = c("Constant", "post"),
    dep.var.labels = dvs,
    covariate.labels = c("(Intercept)", "Post-invasion period"),
    omit.stat = c("f", "ser", "n"),
    column.sep.width = "0pt",
    dep.var.caption = "",
    font.size = "small",
    digits = 3,
    digits.extra = 0,
    star.char = c("+", "*", "**", "***"),
    star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
    notes.align = "c",
    align = TRUE,
    no.space = TRUE,
    notes = "\\parbox[t]{0.9\\textwidth}{\\footnotesize \\textit{Note:} $+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001$. Observations in control group weighted by propensity score, using optimal full matching on age, income, ideology, gender, education, and social class.}",
    notes.label = "",
    add.lines = list(
      c("Indiv-level covariates", rep("No", 5)),
      c("Region (CCAA) FE", rep("No", 5)),
      c("Matched data", rep("Yes", 5)),
      c("Observations", sapply(mlist, function(x) nobs(x))),
      c("Control", sapply(mlist, function(x) extract(x, "c"))),
      c("Treated", sapply(mlist, function(x) extract(x, "t")))),
    notes.append = FALSE),
  filecon)
close(filecon)

# Table (full, with covariates)
filecon = file("tables/tab_barometer_full_covs.tex")
writeLines(
  stargazer(mlist_covs,
    keep = c("Constant", "post"),
    label = "tab:lm_baro_covs",
    title = "Effect of invasion on electoral behavior (no matching, with covariates)",
    order = c("Constant", "post"),
    dep.var.labels = dvs,
    covariate.labels = c("(Intercept)", "Post-invasion period"),
    omit.stat = c("f", "ser", "n"),
    column.sep.width = "0pt",
    dep.var.caption = "",
    font.size = "small",
    digits = 3,
    digits.extra = 0,
    star.char = c("+", "*", "**", "***"),
    star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
    notes.align = "c",
    align = TRUE,
    no.space = TRUE,
    notes = "\\parbox[t]{0.9\\textwidth}{\\footnotesize \\textit{Note:} $+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001$. Individual-level controls include age, income, ideology, gender, education, and social class. Region (CCAA) fixed effects included.}",
    notes.label = "",
    add.lines = list(
      c("Indiv-level covariates", rep("Yes", 5)),
      c("Region (CCAA) FE", rep("Yes", 5)),
      c("Matched data", rep("No", 5)),
      c("Observations", sapply(mlist_covs, function(x) nobs(x))),
      c("Control", sapply(mlist_covs, function(x) extract(x, "c"))),
      c("Treated", sapply(mlist_covs, function(x) extract(x, "t")))),
    notes.append = FALSE),
  filecon)
close(filecon)

# Table (excluding February 24)
filecon = file("tables/tab_barometer_no24.tex")
writeLines(
  stargazer(mlist_no24,
    se = starprep(list(m1b, m2b, m3b, m4b, m5b),
      clusters = matched_data_no24$subclass, se_type = "stata"),
    keep = c("Constant", "post"),
    label = "tab:lm_baro_no24",
    title = "Effect of invasion on electoral behavior (excluding February 24th)",
    order = c("Constant", "post"),
    dep.var.labels = dvs,
    covariate.labels = c("(Intercept)", "Post-invasion period"),
    omit.stat = c("f", "ser", "n"),
    column.sep.width = "0pt",
    dep.var.caption = "",
    font.size = "small",
    digits = 3,
    digits.extra = 0,
    star.char = c("+", "*", "**", "***"),
    star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
    notes.align = "c",
    align = TRUE,
    no.space = TRUE,
    notes = "\\parbox[t]{0.9\\textwidth}{\\footnotesize \\textit{Note:} $+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001$. Treated sample excludes observations from February 24th. Observations in control group weighted by propensity score, using optimal full matching on age, income, ideology, gender, education, and social class.}",
    notes.label = "",
    add.lines = list(
      c("Indiv-level covariates", rep("No", 5)),
      c("Region (CCAA) FE", rep("No", 5)),
      c("Matched data", rep("Yes", 5)),
      c("Observations", sapply(mlist_no24, function(x) nobs(x))),
      c("Control", sapply(mlist_no24, function(x) extract(x, "c"))),
      c("Treated", sapply(mlist_no24, function(x) extract(x, "t")))),
    notes.append = FALSE),
  filecon)
close(filecon)

# Table (excluding February 24)
filecon = file("tables/tab_barometer_no24_covs.tex")
writeLines(
  stargazer(mlist_no24_covs,
    keep = c("Constant", "post"),
    label = "tab:lm_baro_no24_covs",
    title = "Effect of invasion on electoral behavior (excluding February 24th, no matching, including covariates)",
    order = c("Constant", "post"),
    dep.var.labels = dvs,
    covariate.labels = c("(Intercept)", "Post-invasion period"),
    omit.stat = c("f", "ser", "n"),
    column.sep.width = "0pt",
    dep.var.caption = "",
    font.size = "small",
    digits = 3,
    digits.extra = 0,
    star.char = c("+", "*", "**", "***"),
    star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
    notes.align = "c",
    align = TRUE,
    no.space = TRUE,
    notes = "\\parbox[t]{0.9\\textwidth}{\\footnotesize \\textit{Note:} $+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001$. Individual-level controls include age, income, ideology, gender, education, and social class. Region (CCAA) fixed effects included.}",
    notes.label = "",
    add.lines = list(
      c("Indiv-level covariates", rep("Yes", 5)),
      c("Region (CCAA) FE", rep("Yes", 5)),
      c("Matched data", rep("No", 5)),
      c("Observations", sapply(mlist_no24_covs, function(x) nobs(x))),
      c("Control", sapply(mlist_no24_covs, function(x) extract(x, "c"))),
      c("Treated", sapply(mlist_no24_covs, function(x) extract(x, "t")))),
    notes.append = FALSE),
  filecon)
close(filecon)



## Placebo tests and different bandwidths

## Placebo (25-26 vs 27-28)
d_plac = subset(d, fecha %in% 25:28) %>%
  mutate(post = ifelse(fecha > 26, 1, 0))
match_placebo = matchit(f, data = d_plac, method = "full")
matched_data_placebo = match.data(match_placebo)

m_pla1 = lm_robust(pr_voto ~ post,
  data = matched_data_placebo, weights = weights, clusters = subclass, se_type = "stata")
m_pla2 = lm_robust(pr_voto_psoe ~ post,
  data = matched_data_placebo, weights = weights, clusters = subclass, se_type = "stata")
m_pla3 = lm_robust(pr_voto_pp ~ post,
  data = matched_data_placebo, weights = weights, clusters = subclass, se_type = "stata")
m_pla4 = lm_robust(pr_voto_up ~ post,
  data = matched_data_placebo, weights = weights, clusters = subclass, se_type = "stata")
m_pla5 = lm_robust(pr_voto_vox ~ post,
  data = matched_data_placebo, weights = weights, clusters = subclass, se_type = "stata")

m_coefs_plac = bind_rows(
    tidy(m_pla1, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m_pla2, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m_pla3, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m_pla4, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m_pla5, conf.int = TRUE) %>% filter(term == "post")) %>%
  mutate(outcome2 = recode(outcome,
    "pr_voto" = "Would vote",
    "pr_voto_psoe" = "Vote for PSOE",
    "pr_voto_pp" = "Vote for PP",
    "pr_voto_up" = "Vote for UP",
    "pr_voto_vox" = "Vote for VOX")) %>%
  mutate(
    conf.high90 = estimate + std.error * qnorm(0.95),
    conf.low90 = estimate - std.error * qnorm(0.95))

main = ggplot(m_coefs_plac, aes(x = estimate, y = outcome2)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
    fatten = 1, size = 2, shape = 21, fill = "white") +
  geom_pointrange(aes(xmin = conf.low90, xmax = conf.high90),
    fatten = 1, size = 3, shape = 21, fill = "white") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  labs(x = "\nCoefficient estimate of placebo date", y = "",
    subtitle = "Bars indicate 90% and 95% CIs.") +
  theme(axis.text.y = element_text(size = 12))
ggsave("figures/barometer_placebo.pdf", width = 5, height = 3, device = cairo_pdf)

## Compare all days
d_2223 = subset(d, fecha %in% 22:23) %>% mutate(post = ifelse(fecha == max(fecha), 1, 0))
match_2223 = matchit(f, data = d_2223, method = "full")
matched_data_2223 = match.data(match_2223)
d_2324 = subset(d, fecha %in% 23:24) %>% mutate(post = ifelse(fecha == max(fecha), 1, 0))
match_2324 = matchit(f, data = d_2324, method = "full")
matched_data_2324 = match.data(match_2324)
d_2425 = subset(d, fecha %in% 24:25) %>% mutate(post = ifelse(fecha == max(fecha), 1, 0))
match_2425 = matchit(f, data = d_2425, method = "full")
matched_data_2425 = match.data(match_2425)
d_2526 = subset(d, fecha %in% 25:26) %>% mutate(post = ifelse(fecha == max(fecha), 1, 0))
match_2526 = matchit(f, data = d_2526, method = "full")
matched_data_2526 = match.data(match_2526)
d_2627 = subset(d, fecha %in% 26:27) %>% mutate(post = ifelse(fecha == max(fecha), 1, 0))
match_2627 = matchit(f, data = d_2627, method = "full")
matched_data_2627 = match.data(match_2627)
d_2728 = subset(d, fecha %in% 27:28) %>% mutate(post = ifelse(fecha == max(fecha), 1, 0))
match_2728 = matchit(f, data = d_2728, method = "full")
matched_data_2728 = match.data(match_2728)

m1_2223 = lm_robust(pr_voto ~ post, data = matched_data_2223,
  weights = weights, clusters = subclass, se_type = "stata")
m1_2324 = lm_robust(pr_voto ~ post, data = matched_data_2324,
  weights = weights, clusters = subclass, se_type = "stata")
m1_2425 = lm_robust(pr_voto ~ post, data = matched_data_2425,
  weights = weights, clusters = subclass, se_type = "stata")
m1_2526 = lm_robust(pr_voto ~ post, data = matched_data_2526,
  weights = weights, clusters = subclass, se_type = "stata")
m1_2627 = lm_robust(pr_voto ~ post, data = matched_data_2627,
  weights = weights, clusters = subclass, se_type = "stata")
m1_2728 = lm_robust(pr_voto ~ post, data = matched_data_2728,
  weights = weights, clusters = subclass, se_type = "stata")

m2_2223 = lm_robust(pr_voto_psoe ~ post, data = matched_data_2223,
  weights = weights, clusters = subclass, se_type = "stata")
m2_2324 = lm_robust(pr_voto_psoe ~ post, data = matched_data_2324,
  weights = weights, clusters = subclass, se_type = "stata")
m2_2425 = lm_robust(pr_voto_psoe ~ post, data = matched_data_2425,
  weights = weights, clusters = subclass, se_type = "stata")
m2_2526 = lm_robust(pr_voto_psoe ~ post, data = matched_data_2526,
  weights = weights, clusters = subclass, se_type = "stata")
m2_2627 = lm_robust(pr_voto_psoe ~ post, data = matched_data_2627,
  weights = weights, clusters = subclass, se_type = "stata")
m2_2728 = lm_robust(pr_voto_psoe ~ post, data = matched_data_2728,
  weights = weights, clusters = subclass, se_type = "stata")

m3_2223 = lm_robust(pr_voto_pp ~ post, data = matched_data_2223,
  weights = weights, clusters = subclass, se_type = "stata")
m3_2324 = lm_robust(pr_voto_pp ~ post, data = matched_data_2324,
  weights = weights, clusters = subclass, se_type = "stata")
m3_2425 = lm_robust(pr_voto_pp ~ post, data = matched_data_2425,
  weights = weights, clusters = subclass, se_type = "stata")
m3_2526 = lm_robust(pr_voto_pp ~ post, data = matched_data_2526,
  weights = weights, clusters = subclass, se_type = "stata")
m3_2627 = lm_robust(pr_voto_pp ~ post, data = matched_data_2627,
  weights = weights, clusters = subclass, se_type = "stata")
m3_2728 = lm_robust(pr_voto_pp ~ post, data = matched_data_2728,
  weights = weights, clusters = subclass, se_type = "stata")

m4_2223 = lm_robust(pr_voto_up ~ post, data = matched_data_2223,
  weights = weights, clusters = subclass, se_type = "stata")
m4_2324 = lm_robust(pr_voto_up ~ post, data = matched_data_2324,
  weights = weights, clusters = subclass, se_type = "stata")
m4_2425 = lm_robust(pr_voto_up ~ post, data = matched_data_2425,
  weights = weights, clusters = subclass, se_type = "stata")
m4_2526 = lm_robust(pr_voto_up ~ post, data = matched_data_2526,
  weights = weights, clusters = subclass, se_type = "stata")
m4_2627 = lm_robust(pr_voto_up ~ post, data = matched_data_2627,
  weights = weights, clusters = subclass, se_type = "stata")
m4_2728 = lm_robust(pr_voto_up ~ post, data = matched_data_2728,
  weights = weights, clusters = subclass, se_type = "stata")

m5_2223 = lm_robust(pr_voto_vox ~ post, data = matched_data_2223,
  weights = weights, clusters = subclass, se_type = "stata")
m5_2324 = lm_robust(pr_voto_vox ~ post, data = matched_data_2324,
  weights = weights, clusters = subclass, se_type = "stata")
m5_2425 = lm_robust(pr_voto_vox ~ post, data = matched_data_2425,
  weights = weights, clusters = subclass, se_type = "stata")
m5_2526 = lm_robust(pr_voto_vox ~ post, data = matched_data_2526,
  weights = weights, clusters = subclass, se_type = "stata")
m5_2627 = lm_robust(pr_voto_vox ~ post, data = matched_data_2627,
  weights = weights, clusters = subclass, se_type = "stata")
m5_2728 = lm_robust(pr_voto_vox ~ post, data = matched_data_2728,
  weights = weights, clusters = subclass, se_type = "stata")

m_coefs_days = bind_rows(
    tidy(m1_2223, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m1_2324, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m1_2425, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m1_2526, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m1_2627, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m1_2728, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2_2223, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2_2324, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2_2425, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2_2526, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2_2627, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m2_2728, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3_2223, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3_2324, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3_2425, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3_2526, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3_2627, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m3_2728, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4_2223, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4_2324, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4_2425, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4_2526, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4_2627, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m4_2728, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5_2223, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5_2324, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5_2425, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5_2526, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5_2627, conf.int = TRUE) %>% filter(term == "post"),
    tidy(m5_2728, conf.int = TRUE) %>% filter(term == "post")) %>%
  mutate(outcome2 = recode(outcome,
    "pr_voto" = "Would vote",
    "pr_voto_psoe" = "Vote for PSOE",
    "pr_voto_pp" = "Vote for PP",
    "pr_voto_up" = "Vote for UP",
    "pr_voto_vox" = "Vote for VOX")) %>%
  mutate(model = rep(c("Feb 23 vs Feb 22", "Feb 24 vs Feb 23",
      "Feb 25 vs Feb 24", "Feb 26 vs Feb 25",
      "Feb 27 vs Feb 26", "Feb 28 vs Feb 27"), 5)) %>%
  mutate(modelf = factor(model)) %>%
  mutate(
    conf.high90 = estimate + std.error * qnorm(0.95),
    conf.low90 = estimate - std.error * qnorm(0.95))

days = ggplot(m_coefs_days, aes(x = estimate, y = reorder(model, desc(model)))) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
    fatten = 1, size = 2, shape = 21, fill = "white") +
  geom_pointrange(aes(xmin = conf.low90, xmax = conf.high90),
    fatten = 1, size = 3, shape = 21, fill = "white") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  labs(x = "\nCoefficient estimate of later day", y = "",
    subtitle = "Bars indicate 90% and 95% CIs.") +
  theme(axis.text.y = element_text(size = 12)) +
  facet_wrap(~reorder(outcome2, desc(outcome2)), ncol = 5)
ggsave("figures/barometer_effect_days.pdf", width = 10, height = 3, device = cairo_pdf)
